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Abstract: 

This is a paper about multi-fractal scaling and dissipation in a shell model of tur- 
bulence, called the Gledzer-Ohkitani-Yamada model or GOY model. This set of 
equations describes a one dimensional cascade of energy towards higher wave vec- 
tors. When the model is chaotic, the high-wave- vector velocity is a product of 
roughly independent multipliers, one for each logarithmic momentum shell. The 
appropriate tool for studying the multifractal properties of this model is shown to 
be the energy flux on each shell rather than the velocity on each shell. Using this 
quantity, one can obtain better measurements of the deviations from Kolmogorov 
scaling (in the GOY dynamics) than were available up to now. These deviations are 
seen to depend upon the details of inertial-range structure of the model and hence 
are not universal. However, once the conserved quantities of the model are fixed 
to have the same scaling structure as energy and helicity, these deviations seem to 
depend only weakly upon the scale parameter of the model. We analyze the connec- 
tion between multifractality in the velocity distribution and multifractality in the 
dissipation. Our arguments suggest that the connection is universal for models of 
this character, but the model has a different behavior from that of real turbulence. 
We also predict the scaling behavior of time correlations of shell- velocities, of the 
dissipation, and of Lyapunov indices. These scaling arguments can be carried over, 
with little change, to multifractal models of real turbulence. 

PACS: 47.27.-i, 47.27.Jv, 05.45. +b 
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1 Introduction 



The recent literature contains two alternative views of the nature of well-developed 
turbulence. In one view, the simple scaling caught by the K41 [1] paper is the 
asymptotic truth which holds in the limit of high Reynolds numbers. Then the 
experimental facts, some of which seems to support a more complicated scaling, are 
described in terms of nonasymptotic corrections to scaling [2, 3, 4, 5, 6, 7, 8, 9]. In 
the other view [10] the experiments are better understood and described as a result 
of a multifractal picture [11, 12, 13] in which cascades produce anomalously large 
fluctuations in the velocity fields. These two views can both be supported by the 
experimental evidence [14, 15, 16, 17, 18, 19]. There are theoretical arguments for 
both. 

It would be ideal to distinguish these two possibilities by direct numerical sim- 
ulations of Navier- Stokes dynamics. Unfortunately, the current computing power 
sets an upper limit on the Reynolds number. So far, the highest Taylor Reynolds 
number is on the order of 200 [20, 21], which is insufficient to make the distinction. 
Consequently, many people have worked on simplified approaches which might offer 
some understanding of Navier-Stokes dynamics. 

One approach, called the reduced wave vector set approximations (REWA), ap- 
proximates the Navier-Stokes dynamics by representing the full velocity field on a set 
of wave vectors which gets more and more thinned out for higher wave vectors. For 
detailed discussions of the method and the results we refer to refs. [22, 23, 24, 2, 3]. 

Earlier, Gledzer [25] introduced another, and simpler representation of Navier- 
Stokes dynamics in which the velocities are placed on a one dimensional array of 
wave vectors. Each successive new velocity falls in a new shell in k-space in which 
the wave vector is increased by a factor of A so that the nth shell has k n = k X n . In 
the version of the model [26, 27, 28, 29, 30] used here and referred to as the GOY 
model, each shell is described by a single complex variable, the complex velocity U n . 

The GOY model shows at least two qualitatively different kinds of behavior [31]. 
In one range of parameters, the system relaxes to a time-independent state in which 
the velocity decays (apart from boundary corrections) with wave vector according 
to the predictions of K41. However, for other parameter values, the system has a 
long-term behavior which includes stochastic fluctuations in the velocity. The basic 
statistical variable is the ratio of the velocity fluctuations in neighboring shells. If 
this multiplier fluctuates locally, there cannot be any locking of the correlations 
among the fluctuations in far-distant shells. There is no way that independent 
short range finite strength interactions can, in one dimension, be translated into 
infinitely strong long ranged correlations. (The argument for this special behavior 
of one-dimensional systems was originally given by Landau [32] in the context of 
phase transition problems. Of course in higher dimensions short ranged interactions 
can indeed produce long-ranged correlations via phase transitions, see below.) The 
argument about relatively weak correlations in one dimension was first applied to the 
GOY model by Benzi, Biferale and Parisi [30] (hereafter referred to as BBP). They 
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argued, whenever there are any fluctuations in the long time dynamics, the GOY 
model necessarily develops very considerable fluctuations in the ratio of velocities in 
far-distant shells. As a result, the velocities will show a multifractal behavior in the 
inertial range and there will be correlations among far distant shells [33, 34, 28, 30, 
29, 35]. In contrast, REWA has its dynamics on a higher dimensional ^-lattice. The 
dynamics can then possibly develop long range correlations among velocities with 
different A;- vectors. These correlations may damp down the largest fluctuations in 
the velocity amplitudes and produce the observed asymptotic behavior [24, 3] which 
is fractal (K4I) rather than multi-fractal. 

Compared to REWA, the simpler GOY model offers the opportunity of a some- 
how easier analytic approach, and of course of less expensive numerics. REWA in 
turn offers the advantages of additional degrees of freedom and of an apparent extra 
closeness to Navier Stokes. But, we cannot tell which approach comes closer to the 
truth of the full Navier Stokes system. 

In this paper we aim for the more modest goal of finding some of the properties 
of the GOY model. Specifically we wish to 

1. Develop a more accurate numerical method of determining the scaling ex- 
ponents of the multifractal spectrum and thereby to better understand the 
spectrum of fluctuations of \U n \. 

2. Understand the role of conserved quantities. 

3. Examine whether the multifractal properties dependent upon the exact form 
of the dissipation and of the nonlinear term. 

4. Describe and verify a theory for time dependence of various quantities of the 
dissipation and thereby generate a scaling theory for the multifractal properties 
of the dissipation, for velocity correlations, and for fluctuations in Lyapunov 
indices. 

In the next section of this paper, we describe the GOY model and some of its 
major qualitative features. Section three is devoted to the properties of the velocity 
in the inertial range. We show how to get more accurate values of the multifractal 
scaling exponents than were available up to now. We then see how the exponents 
depend upon the various parameters of the model. In section four, we do some 
calculations related to time-dependent correlations in the model: subsection one 
describes the relation between the multifractal properties of the dissipation and 
that of the velocity in the inertial range; the next estimates the order of magnitude 
of velocity correlations; and the last discusses fluctuations in the Lyapunov index 
for the model. 
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2 The GOY model of turbulence 



The GOY model describes a one-dimensional cascade of energies among a set of 
complex velocities, U n , on a one dimensional set of wave vectors 

k n = k X n , n = 1,2, ...,7V. (1) 

The model is a system of ODE with the following structure: 

^-U n = -D n + F n + iC n . (2) 
at 

Here D n stands for a dissipation term, F n stands for a forcing term which is only set 
on low-n shells, and C n stands for nonlinear couplings among different shells. The 
last term is crucial to inducing energy cascades in the model. 

We shall make the same choices for the three terms as have been used in several 
previous turbulence studies[25, 26, 27, 28, 30, 29, 31] : 



D n = vk 2 n U n , (3) 
F n = fS nA , (4) 
C n = aKU^^ + bK^^ + cK-^UZ^. (5) 

They are intended to capture some of the features in hydrodynamics: viscous dis- 
sipation 1 of energy, external forcing on a large scale, and quadratic interaction 
among different modes with strength proportional to k. Furthermore, we impose 
rigid boundary conditions on U n in which the only non-zero U n 7 s are those for which 
n is within the range [l,iV]. The constants, a, 6, c, /, A, and k define the model. 
Throughout this paper we make the conventional choices 



k = 2" 4 , (6a) 

/ = 5(1 + - 10" 3 , (6b) 

a = 1, (6c) 

b = -e, (6d) 

c = -1 + e. (6e) 

The standard case which we will use for comparison with the results in the literature 
has 

A = 2, (7a) 

v = 10" 7 , (7b) 

e = 1/2, (7c) 

N = 22. (7d) 



1 Later on, we will also allow for hyperviscosity D n = vkfjj n and study the behavior at <f> = 4, 6. 
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Such a system is similar to the 3d Navier-Stokes dynamics in four respects: 

1. In the inviscid and non-forcing limit, there are two conserved quantities [25] 
which can be identified with the total energy / \u(x) | 2 dx/2 and the helicity 
f u • V X udx of the exact problem. 

2. The cascade term conserves the phase volume, defined as the total volume in 
the iV-dimensional complex velocity space. The result is a direct consequence 
of the statement that 

m c - = °- (8) 

3. The system can reach a steady state in which it behaves chaotically. Since 
the system is forced at large scale and the dissipation occurs mostly at small 
scales, the system must cascade energy from large scales to smaller ones. 

4. The multifractal behavior shows some resemblance to the behavior seen ex- 
perimentally. 

A general discussion of conserved quantities can be found in Gledzer's work [25]. 
We simply notice that the cascade terms give rises to the expression for the general 
conserved quantity: 

W = J2\U n \ 2 z n } (9) 

n 

whenever z satisfies the quadratic equation 

= a + bz + cz 2 . (10) 

We require one of the two conserved quantities to have the structure of the kinetic 
energy, i.e. z = 1: 

£ = El^l 2 /2. (11) 

n 

Then from equation (10), we must take a + b + c = 0, as reflected in equations (6c) 
- (6e). By adjusting the time scale, we can make a = 1. So the cascade terms of the 
model contain only two free parameters, which can be defined to be e and A. Now 
the energy conservation law takes the form 

±\ Un \2 = -vk 2 n \U n \ 2 + nfU^n,4] + Jn-1 ~ Jn- (12) 

at 

Here the different terms refer respectively to the dissipation, to the forcing, and to 
a discrete divergence of an energy flux. The cascade produces the fluxes, which are 
defined in terms of the triple products: 

A n = k n -iU n -\U n U n +\. (13) 

Then the energy flux from the nth mode to the n + 1th mode is 

J„ = ^[-A„ +1 -(l-e)A„]. (14) 
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We therefore can picture a steady state of the dynamical system as a cascade of 
energy from large "eddies" to smaller ones, where the energy is dissipated through 
viscous diffusion. It is in this sense that we say the dynamics may mimic real 
turbulence. 

Given an e, the other conserved quantity is then determined by equation (9) in 
terms of the other solution of equation (10), z = l/(e — I) as 

H = Y,\U n \ 2 (e-l)- n . (15) 

n 

The corresponding conservation law takes the form of 



4(e - iy n \U n \ 2 = -vkl\U n \\e - l)~ n + K[fU:S n4 ](e - l)~ n + L n . x - L n , (16) 
at 



where 



L n = (e-l)- n Z[(A n -A n+1 )] (17) 



is the relevant flux from the (n — l)th mode to the nth mode. 

The inviscid three dimensional hydrodynamics has two conserved quantities 
quadratic in velocity: the total energy and the helicity; the latter is a spatial integral 
of velocity dotted into vorticity. In our shell model, the closest analog to helicity is 
J2 n (~ ^) n k n \U n \ 2 . Whenever 1/(1 — e) is equal to A, the helicity defined in this way 
agrees with our second conserved quantity H . In particular, this equality holds for 
the conventional choice of parameters (6) with e = 0.5. 

We remark that the second conserved quantity plays an essential role in the 
theories [37, 3, 4, 38] which use corrections to scaling to explain why K4I does not 
fully fit experimental turbulence data. In their approach, the leading correction to 
scaling is the result of an additive correction to the asymptotic scaling. 

The structure of the cascade part of the equation of motion is particularly inter- 
esting. This part has the form 

r -V/ — — (^) 
n ~h su m su p - 

Here E and H are our two conserved quantities and J is proportional to a completely 
antisymmetric function of its three indices. The equation of motion generated from 
equation (18) 

^-U n = {U n ,E} with (19a) 

at 

r , , 6A SB 6H , , , 

{a, b} = i J2 J ™7?r77r77r + cx - ( 19b ) 

n,m,p 0U n dU m dUp 

looks as if it might be Hamiltonian in character in that it is generated by an anti- 
symmetrical bracket structure. The form given by equation (19b) suggests that what 



7 



we are seeing is possibly a Lie-Poisson System [39]. However, detailed investigations 
by one of us (MM) have shown that this bracket in particular and any other possible 
Lie-Poisson-like bracket generating the equations of motion fails to satisfy the Jacobi 
identity, so that the brackets in question are not really Poisson brackets. 

The numerical integration of the shell equation (2) was done using the 'lsode'- 
software package [40] , a differential equation solver that handles the stiffness of the 
equations efficiently. We benchmark our code by running a well studied case of 
equation(7) [29, 28, 26], in which the total number of shells is 22. Our results agree 
with the previous ones. The satisfaction of the balance eqs. (12) and (16) for the 
stationary state provides an extra check for our numerics. 

3 Wiggles in the average of \U n \ 
3.1 Static Solutions 

Biferale et al. [31] studied two classes of solutions to the GOY model. They observed 
numerically (for A = 2 and v = 10~ 7 ) that when < e < 0.3843, the system 
reaches a stable static solution of the Kolmogorov type. On the other hand, when 
0.3953 < e < 2 their analysis shows a chaotic and time-dependent behavior. Much 
of the previous analysis of the GOY model has been done for e = 0.5. Since the 
phase transition at e ~ 0.4 apparently has the character of a continuous transition, 
and since e = 0.5 is reasonably close to 0.4, we should expect many of the qualitative 
features of the static state to manifest themselves for e = 0.5. 

To understand the structure of the nonlinear interaction, Biferale, et a/.[31] ana- 
lyzed the static solutions in the inviscid and homogeneous GOY model to obtain an 
iteration map for the ratios of velocity. Here we offer a slightly different approach 
which gives an analytic solution for the ratios in the inhomogeneous and inviscid 
case. Let us first consider the case in which v = and / = 0. For n > 4, one 
sets the cascade operator of equation (5) equal to zero and finds a linear difference 
equation for the product of three velocities, A n , defined by equation (13): 

A n+1 - eA n - (1 - e)A n _ 1 = 0. (20) 

This equation admits a general solution of the form: 

A n = A + B(e-l) n } (21) 

where A and B are each complex constants, related to the fluxes for energy and the 
other conserved quantity by 

J n = 3?[(e - 2) A] = 3ft(/*7 4 *) (22a) 
L n = Z[(e-2)B} = %(fU:)(e-l)-* (22b) 

Notice that, in the static case, the ratios of U n } s can be expressed in terms of the 
ratios of A n 's, in particular, 
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r n = U n+3 /U n = A n+2 /(XA n+1 ). (23) 

Rewriting equation (20) in terms of r n gives the ratio map obtained by Biferale et. 
al.. The general expression (21) permits two simple scaling solutions corresponding 
to B = and A = 0, respectively. However, if A and B are both of order unity, 
then, as long as < e < 2, the first scaling solution (B = 0) will always dominate 
for large re. 

For the case B = 0, A ^ 0, we have A n = k n _iU n -\U n U n+ \ = A, independent of 
re, and r n = 1/A. Note that the velocity amplitudes show period three oscillations 
superimposed on the classical k~ x l 3 falloff, i.e. 

{Uokn -1 / 3 for re = (mod 3) 
u x k n ~ xl3 for re = 1 (mod 3) (24) 
U2k n ~ 1 / 3 for re = 2 (mod 3) 

The period three oscillations will be important in what follows. 

For the second scaling solution (A = 0, B ^ 0) we have r n = (e — 1)/A. In 
particular for A = 1/(1 — e) (e.g., A = 2, e = 1/2), where L n can be considered as 
helicity flux, this results in \U n \ oc k~ 2 / 3 [37, 3, 4], again, up to a period 3. 

Let us now allow for the inhomogenity in the difference eqs. (20) by including 
the forcing term fU*6 n 4 } but still neglecting viscosity. In the limit as viscosity goes 
to zero, Uz also goes to zero. Consequently, for the time-independent situation, we 
find A 4 = and then, A/B = (e — 1)~ 4 . In this way, we derive 

l-(e-l) n - 2 , , 

rn = x[1 _\ e _; )n -sy - > 3. (25) 

This expression for r n has an oscillatory behavior for re in the inertial range whenever 
e < 1. Note that this solution is determined by our choices of boundary conditions 
and the stirring mechanism. 

In table 1 we show the comparison of the numerical results and the analytical 
expression (25) for r n in the inertial range when e = 0.3 (static solution) and e = 0.5 
(chaotic solution where we define r n = (U n+ 3 )/(U n )) } respectively. In the former 
case, the agreement between the numerics and the theory is very good, with small 
deviations resulting from neglecting the viscous term in (25). In the latter case, 
we notice that although the values do not match closely, nevertheless the static 
theory (from which we derived (25)) describes correctly the oscillatory behavior of 
r n . This suggests that some of the quantities in the static solutions may well reflect 
the properties of chaotic solutions. 

3.2 Velocities, Triple products, and Fluxes 

Previous works [28, 29, 30] studied moments of the U n } s: 

s n , q = {\U n \ q )^k-^. (26) 
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Numerically, we observe that the scaling of the U n } s exhibits, superimposed on the 
overall power law scaling, the period three oscillations which have been reported 
earlier by Pisarenko et. al.[29]. The phenomenon is reminiscent of and probably has 
the same origin (cf. eq. (24)) as the period three oscillations in the static Kolmogorov 
state. These oscillations are shown in Figure 1 for two cases with different forms 
of dissipation: a normal viscosity case (D n = vk 2 n ll n ) vs. a hyperviscosity case 
(D n = vk^Un). Oscillations are present in both cases, but much more pronounced 
in the latter situation. 

We could improve our scaling analysis if we did not have the period three oscil- 
lation which confuses the pictures of power law decay. In the static case a product 
of the form (s^i^^^^i^) 1 / 3 would not show any oscillation. To verify that this 
result extends to the chaotic situation, we plotted (s^i^Sn^s^i^) 1 / 3 . Indeed, the 
oscillations disappear. 

The period three oscillations are responsible for the large uncertainties in mea- 
suring the scaling exponents. This uncertainty can be eliminated by extending our 
analysis of the static state in the previous section. We first observe that relations 
similar to (20) (21) (22b) hold for 5(A„), and we obtain: 

3(A n ) = K^U^UnU^) = -^±L (l - ( £ - l)"" 4 ) , (27) 

We are then motivated to study the scaling of the triple moments 

S n , q = {\ < 3[U n - 1 U n U n+1 ]\ q ' 3 ) oc k~<«, (28) 

which are free of period three oscillations. As one can see from Figure 2 these 
quantities do show a cleaner scaling form than the s n ^ q . 

However, the oscillations at small n still remain a problem. We then make use 
of the analog to the Kolmogorov structure equation [33] 

(J n ) = -%((A n+1 ) + (1 - e)(A„)) = ®(m)- (29) 

In the inertial range, J n should be independent of n. Eq. (29) suggests that we could 
eliminate small-n oscillations by studying the scaling of 



Sn,g — ( 



q/3 

>, (30) 



U n U n+ iU n+ 2 + ^ — - — ^ U n -\U n U n+ \ 
which, in mean field approximation, scales as 

Xn, q OC k^ 3 . (31) 

Expression (31) has the same scaling behavior as S n>q , but is free from the oscillation 
for small n in the inertial range. 

This sequence of operations enable us to reduce the statistical error in 6( q from 
0.05 down to 0.005, and to increase the scaling region by perhaps 1 or 2 decades. 
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The improved accuracy allows us to study the dependence ol scaling exponents on 
various parameters in the model such as e, and the form of dissipation 2 . 

Note that the comparably large spectral strength s n for n = 17 (cf. fig. 1), i.e., 
just before the viscous range sets in, can be understood as a kind of bottleneck 
phenomenon [5, 6]. As pointed out in ref. [6], the bottleneck energy pileup is more 
pronounced for the hyperviscous case, since viscosity sets in more rapidly. On the 
other hand, there cannot be a bottleneck energy pileup for the quantities S n and 
E n , because these quantities are based on the energy flux, which is constant in the 
inertial range. 

3.3 Multiscaling 

We have shown that by analyzing first s„ ig , then S n>q , and finally S„ jg , the accuracy 
of the scaling exponents ( q can be improved step by step 3 . Besides the statistical 
errors reported above, the analytical expressions (27) and (29), which state ( 3 = 1, 
help us to control the systematic error in the ( q determination, since consequently 
S n> 3 and E„ j3 should scale as oc k~ x . The numerical value for the scaling exponent 
of E„ j3 is indeed very close to 1, e.g., ( 3 = 1.005. We choose the averaging time for 
each simulation by demanding that S n>q , and E„ jg have all settled down to their 
long-term value. This stability is ensured by comparing successive runs, for each of 
which the averaging time of 15 000 large eddy turnover times turns out to be more 
than sufficient for N = 22. The inertial range is determined by the condition that ( 3 
is very close to 1. The numerical values for the deviations from Kolmogorov scaling 
6( q = ( q — q/3 presented below are measured from E„ jg . 

In order to compare the 6( q for different parameters and to make claims on what 
they depend it is rather important to perform a very careful analysis of possible 
errors of our determinations. Besides unknown systematic errors there are three 
controllable source of errors: (i) The inertial range (and thus the fitting range for 
the straight line fit) is not unambiguously determined by the condition "(3 very close 
to 1". In practice, ( 3 is typically off from 1 by an amount of order ±0.005 for typical 
fits in the n-ranges from 4 to 15 or 5 to 15. We attempt to minimize this source 
of error by dividing all measured ( q by ( 3 . Nevertheless, the ( q slightly depend on 
the fitting range, (ii) There is some statistical error from the least square fit of the 
data to a straight line, (iii) Results may differ from run to run, though each is over 
15 000 to 50 000 large eddy turnover times. 

The errors from (i) and (ii) are about of the same size, whereas that from (iii) 
is clearly smaller. So, with error propagation, the total standard deviation is about 
1.5 times as large as the statistical error (ii). These are the errors given in figs. 3-5. 

To visualize the differences in the 6( q for different sets of parameters, we found it 

2 Note that Biferale et al. [31] eliminated the period 2 oscillations by changing the boundary 
conditions of (2). 

3 After completing much of this work, we were told by Z.S. She that his group had independently 
derived a similar method for analysis of the multifractal properties of the GOY model. 
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convenient to plot S( q /(q(q — 3)) rather than 6( q itself. The reason is that we want 
to eliminate the trivial agreement of the 6( q for q = 3 and q = 0, i.e., 6( 3 = 6( = 
for all kinds of parameter sets. 

Figure 3 shows the 6( q and their errors - calculated and displayed as described 
above - for the standard GOY model parameters, A = 2, e = 1/2. In the same figure 
we also give the 6( q obtained from She and Leveque's (SL) [41] result 



( q = 9/9 + 2 



2\ 9/3 

i : 

3 



(32) 



which is picked to be a simple phenomenological description of a situation in which 
one has an upper cutoff to the size of the turbulent fluctuations. Note that in any 
(finite time) numerical calculation such an upper bound will be given. Their result is 
in good agreement with experiment. For q > 3 we find surprisingly good agreement 
also with the numerical GOY results at the standard parameter value. This means 
that the tails of the velocity probability distribution function (PDF), corresponding 
to large fluctuations, are well described by the SL theory. For q < 3 we find slightly 
worse agreement, i.e., the small fluctuations (peak of PDF) are not so well described 
by SL. 

However, the GOY-model values of 6( q are not universal. They depend on the 
choice of e and A. For example, for e = 0.3 and A = 2, 6( q shows classical Komolgorov 
scaling, see Biferale et. al.[31]. Our numerical analysis shows that even when we 
are in a chaotic situation, 6( q varies with A and e. As seen in figure 3, e.g. for 
e = 0.7, A = 2 the scaling corrections are much larger than for the standard case. So 
why should the parameter values A = 2,e = 0.5beso special as to agree with She 
and Leveque? We recall our remark made in section 2 about the second conserved 
quantity. As we have pointed out, given the above special parameters, the extra 
conserved quantity closely resembles helicity, which is indeed a conserved quantity in 
3D inviscid hydrodynamics. We also believe that specific forms of conservation laws 
will have crucial influence on the outcome of statistical averages of the dynamical 
quantities. Consequently, we hypothesize that the scaling corrections should be 
invariant along a curve in e — A plane, which conserves both energy and helicity. 
The curve defined by the above constraints takes the form 

A = r l-. (33) 

To test this hypothesis, we experiment numerically with systems of different (e, A). 
In Figure 4 we present the scaling corrections 6( q for three pairs of e and A 'on the 
curve' (33). And as a comparison, we show a case for which the parameter- values 
are 'off the curve': (e = 0.7, A = 2). The cases 'on the curve' are grouped closely 
together, while the case 'off the curve' is much further away. Our numerical findings 
suggest that the scale factor A plays a relatively minor role in the dynamics once the 
scaling for the second conserved quantity is picked (via equation(33)) to correspond 
to scaling for the helicity. Returning to our original question, we are ready to claim 
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that the canonical choice ol (e = 0.5, A = 2) is a special one because it lies on 
the energy-helicity curve (i.e., respects both the energy and helicity as conserved 
quantities), but it is not a unique one because there are many other parameters on 
the same curve which will give roughly the same values of 6( q . 

After having compared that curve with the She-Leveque theory [41], it is also 
interesting to compare it with the approximations described in BBP. There are two 
steps of approximation in the that paper. In the first, described in section 3 of BBP, 
an approximate recursion relation relating different 6( q is derived. These recursions 
can be solved if one knows the results at q = 1 and q = 2. Using these values as 
adjustable parameters, we have derived all the integer-index 6( q . The parameters 
can be adjusted so that these results agree within our numerical errors with the 
GOY-model numerical solutions. Moveover, in section 4 of BBP, their theory is 
extended to give an approximate value of all 6( q from first principles. When we 
check this part of their theory, we find that there are severe differences between 
their solutions and the numerical calculations. We thus conclude that while the 
section 3 results might be accurate, the section 4 results disagree with the facts. 

As suggested in the conclusion of BBP, correlations in the multiplicative pro- 
cess among the shells may play an important role in order to obtain an accurate 
theoretical estimates of the scaling properties using the closure equations. 

3.4 Dependence upon the form of dissipation 

A further interesting issue is whether the scaling corrections 6( q - besides depending 
on A and e - also depend on the magnitude and form of the the viscous damping 
term. In particular, we take the damping term to have the form D n = vk n ^U n 
and then focus upon the issue of whether the scaling results depend upon v and 
<f>. The analysis of BBP [30] suggests that 6( q is determined only by the form of 
the cascade term, since it is determined by short-ranged correlations between the 
different shells. She [42] suggested an alternative viewpoint: Energy fluctuations are 
produced at all length scales and tend to cascade downward toward smaller lengths 
or higher n-values. When they enter the viscous subrange, they see a changed envi- 
ronment, because the viscosity term then effectively enters the equation of motion. 
Depending upon the value and form of the viscosity term, more or less energy might 
be reflected toward lower n- values. Thus the energy would flow through partially a 
direct and partially an inverse cascade. The amount of reflection would determine 
the corrections to scaling indices. The details of the reflection will be determined in 
part by the form and magnitude of viscous damping. 

Let us look at the facts once again by plotting S( q /(q(q — 3)) as a function of 
q for various cases. Now we fix A and e to their standard values and vary v and 
<f>. We determine the 6( q and their errors exactly as pointed out in the previous 
subsection. The various results for the deviations 6( q - though slightly different - 
are reasonably close together, cf. figure 5. The results for different types of viscosity 
differ by about two standard statistical errors. 
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We have also done a preliminary analysis of the //-dependence. In the hyper- 
viscous case the results seem to be rather somewhat dependent on v . For 5 = 6, 
we found, for different values of z/, both stronger and weaker inertial range scaling 
corrections than for the standard s = 2 situation. 

Our numerical analyses seems to hint that 6( q might depend on the type of 
viscosity. However, our result is certainly not accurate enough to be definitive. 
Many previous workers have gotten confused by crossover effects or by corrections 
to scaling. However, it is to be expected that future workers will obtain improved 
results, perhaps sufficient to tie down the errors and settle the issue. 

4 Time correlations 

4.1 Connection between multifractal behavior of U n and 
dissipation 

In this section, we shall move back and forth between the shell model and the 
theory of real turbulence. In particular, we shall wish to contrast the two theories 
of multifractality, which bear many resemblances. However, as we shall see there is 
also a crucial difference, which is rooted in the nature of the two situations. In real 
turbulence, if there is a large scale velocity U 0} then any small scale disturbance 
will produce a time derivative of the velocity 

MM = -LVVL7(JM). (34) 

This time derivative produced by spatial variation and carried by the large-scale 
motion is called sweeping. There is no analog of sweeping in the GOY model since 
there is no coupling between Uq and the large- n U n (t). This distinction will produce 
a considerable difference in the answer which we will obtain for the multifractal 
properties of dissipation. 

Recall that in the multifractal approach to the usual turbulence theory [11], 
averages of the velocity difference at a distance r have the form 

< [U(R + r,t)-U(R,t)] q >~V q [r/L] c \ (35) 

where Vo is a typical value of the velocity, L is an integral scale and ( q is the 
multifractal scaling exponent for the velocity. The < • • • > in equation (35) stands 
for an average over space and time. The analogous quantity within the shell model 
is s njq (or S n>q and E„ jg correspondingly), which scales with wave vector k n , 

s n , q ~ k n ~ Cq , (36) 

where now the average implied is a time average. In both cases, the deviations of ( q 
from q/3 measure the deviation from the K41 description. 
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Next we look at another form ol multilractal behavior. As the energy cascades 
towards higher values of n it reaches sufficiently high wave vectors so that the 
viscosity can produce energy dissipation. According to equation (12) the rate of 
energy dissipation is 4 

e(t) = uY J k n 2 \U n (t)\ 2 . (37) 

n 

An analogous dissipation occurs in the Navier Stokes theory. There the dissipation 
depends upon both space and time and takes the form 

e(R, t) = l -v Y^(d iUj (R) + d jUl (R)) 2 - (38) 

To get an appropriate fluctuating quantity we take the dissipation and average it 
over a time interval r. For the GOY case take the average to be 

e T (t) = J* +T e(t')dt' . (39) 

In contrast, Meneveau and Sreenivasan [43] looked at real turbulence data obtained 
by taking a particular point in space and averaging over a period of time. They 
measured time averages as 

e T (t) = J* +T e(R,t')dt'. (40) 

To be more precise, e(R,t) was substituted by its one dimensional surrogate e'it) ~ 
(d t Ui) 2 . The dissipation is thus measured always at the same position and there is 
no reason to include its independence. In both situations, the multifractal behavior 
of the dissipation is defined by considering averages of powers of the dissipation in 
the form 

< l^W >~ T^. (41) 

If fi q is not proportional to q for r in the inertial range, then the dissipation is said 
to be multifractal. 

Following [10] [11] [48], Benzi et al [36] developed a theory of this multifractal 
dissipation for the case of Navier Stokes turbulence, noticing that the dissipation of 
energy was fed through a cascade in which the flux goes through all scales of r up 
to and including the dissipation scale. 

The scales in space and time are connected by the sweeping process. If the inertial 
scale velocity is of order U then the scales are connected by Taylor's hypothesis [45] 

r ~ (42) 



4 The reader should not mix up the energy dissipation e(t) with the GOY model parameter e 
introduced in eq. (6d). We shall not refer to the latter in this section of the paper. 
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According to the analysis ol [44, 10, 43, 11], \i q is completely determined by ( q . This 
analysis says that to an order of magnitude, the dissipation on scale r is set by the 
flux at that scale, which can be estimated as 

J r (t) ~ r-W^t), (43) 

where U r is the velocity difference on scale r. Putting together equations(41),(42), 
and (43) one finds that the dissipation on scale r has the order of magnitude 

(4) ~< r- q [U r {t)f q >~ T- qH3 ", (44) 

concluding thus that fi q is determined by the multifractal scaling of velocity accord- 
ing to 

Vq = -q + (3q- (45) 

We would like to apply this approach to the shell model. However, this approach 
will not quite work in this case, because there is no sweeping, so Taylor's [45] frozen 
flow hypothesis is not meaningful. Moreover, there is no spatial dependence so one 
can only deal with time averages. Thus in order to perform a similar analysis on 
the shell model, we start again using equation (37) to define the time-integrated 
dissipation. The next ingredient is to get the connection between the shell number, 
n and the time scale r. To an order of magnitude, in the shell model, the time 
derivative of U n is given as k n U n 2 . Thus, there is a natural turnover time for each 
shell defined by 

Tn(t) = IKUnit)]- 1 . (46) 

This equation can be solved to get the dependence of n upon r at each value of t 
and defines the function n(r, t). The final step is to get the dissipation from the 
energy flux. If the shell n is correctly picked so that the time-scale lies within the 
nth shell, then the dissipation should be estimated as the energy flux in that shell. 
In the GOY model then, the dissipation and the flux may be estimated as 

e T (t) ~ Jn(t) ~ ZiKUn-^UnWUn+^t)], (47) 

where in this equation, n is to be evaluated at n(r, t). 

Because U n (t) is a fluctuating quantity, so is T n (t) and we cannot once and for 
all determine which shell belongs to which value of r. This makes the analysis of 
the shell model more difficult than the calculation for the real turbulence problem. 
To analyze the model, we must go back to the f(a) formalism which grew up as 
an alternative description of multifractal phenomena [11, 12, 46]. Here, f(a) is the 
singular spectrum as commonly defined in [12, 46, 43]. Instead of defining ( q , we 
define the probability that U n (t) will have a magnitude which is of the order of 
k n ~ a . This probability is a strongly varying function of n, and of the form k n ^ a \ 
independent of n. Thus, using the delta function one writes 

(6(lnU n (t) + alnk n ))~k f J a ) (48) 
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and calculates averages as integrals over a. Thus for example, 



1=<1>= / da(6(lnU n (t) + alnk n )) ~ / dak n i{a) . (49) 



The integral is to be done by a steepest descent technique. (We only keep factors 
which are exponential function of n. Powers of In k n are neglected.) In this way, 
we learn that the maximum value of f(a) is zero. An analogous calculation gives a 
direct evaluation of ( q . Consider 

< \Un(t)\ q >= [ da(\U n (t)\ q S(lnU n (t) + alnk n )) ~ / dak n ~ qa+s ^. (50) 



On the one hand this average is set by the definition of equation (26) to be k n ~^ q . 
On the other, the integral can be performed via steepest descents and gives the 
familiar Legendre transform definition of f(a) } 

( q = mm[qa - f(a)]. (51) 

Once we know ( q , f(a) is known. 

Next we estimate the average of e T taken to various powers. Once again we 
compute the averages by integrating over a. Now we have two conditions: the first 
being that e T is given by expression (47) for some appropriate value of n, the second 
being that the appropriate value of n is defined via a solution of equation (46). 
Employing these two conditions, we find that 

< e T q >= j da j dn(\J n (t)\ q 6(lnU n {t) + alnk n )6(lnT + lnU n {t)k n )). (52) 

Substituting for the flux using equation (47) and observing that the delta function 
permits us to do the n-integration, we thereby reduce the result to a simple integral 
over a of the form 

< e T q >~ j daT X{a \ (53) 
with the exponent having the value 

X(a) = 3ga " 1 - /(a) . (54) 
1 — a 

In the usual way, the integral is calculated by steepest descent and fi q is determined 
from the saddle point integration as 

fi q = minX(a). (55) 

As desired, this equation connects the scaling exponents fi q and ( q . Note in par- 
ticular, that //i = holds as it should be, independently of the form of f(a). For 
similar approaches, leading to slightly different results, we refer to refs. [48]. 
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The next step is a comparison of our theory with the results from the simulation 
of the model. For convenience, we denote < e q T > by e Tjg . We expect e Tjg to scale 
in the inertial range. It is not clear a priori, what the inertial range will be in 
the time domain, as an application of Taylor's hypothesis in the GOY model is not 
meaningful, see above. We find good scaling behavior for r = 0.2 to r = 6, see Figure 
5 for that of e Tj2 . To extract even preciser scaling exponents from the numerics, we 
plot e Tjg vs. e Tj2 and determine the ratio / u g / / u 2 , which is a direct analogy of extended 
self similarity introduced by Benzi et. al. [47]. Figure 5 shows e Tj4 vs. e Tj2 . Scaling 
is seen in the range between r ~ 0.05 and r ~ 5, where the times are given in the 
natural time units of the GOY model, the large eddy turnover times. 

As stated above, to connect the ( q with the /i g , we need to find the singular 
spectrum f(a) from the ( q . We take advantage of the fact that the numerical values 
coincide with the She-Leveque formula (32) [41]. Then f(a) can be easily obtained 
through the Legendre transformation of ( q , 

F , n 3(a - 1/9) / 

/(a) = -2 + ^ , , V n 
Jy ' In (2/3) V 

We then use our formulas (53), (54) and (55) to find the ratios fi q /fi2- The 
comparison with the numerical values extracted from the scaling of e Tjg versus e Tj2 
is shown in table 2. We estimate the error in the numerical values by comparing 
the results of linear fits in slightly different regions of the scaling range and find it 
to be between 1% for small q and 5% for larger q. Up to the 5 th moments of e T 
the agreement between the numerical results and our scaling analysis seems to be 
within the statistical error. 



1/9 - a 
ln(2/3) 



;i+ln(2/3)) 



(56) 



4.2 Time correlations for \Un 



In this section, we shall estimate the scaling behavior of one of the most basic time 
correlation functions, namely, 



Cp 1 n 1 ,p 2 n2 ( T ) 



{\U ni (t)\ Vl \U n2 (t + r)r) . 



(57) 



We restrict ourselves to the case where r is much larger than the typical time scales 
of the shells n\ and n 2 . The main idea of the calculation is to introduce a (large 
scale) shell rait) into (57) such, that it has a relaxation time of order r, as one 
can then assume, that the shells n\ and n 2 are completely uncorrelated with m, 
whereas there should be more or less full correlation on larger scales than m. We 
thus determine the shell number m through r ~ r m (t), or, by employing (46), 

,(0 = TTT^TH- ~ (58) 



T ~ T„ 



\U m (t)\k r 

with |/7 m (t)| oc k^f^\ We may thus write 

U ni (t) Pl U n2 (t + r) 



U m {t) 



U m {t + T) 



'92 



\u m (t)r\u m (t + t)\ 



V2 
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dm 



u ni (t) 

U m (t) 



pi 



Un 2 (t + T) 
U m (t + T) 



P2 



\u m (t)r\u m (t + t)\ 



P2 



S[(l- f3(t))lnk m + lnr]. 
In terms of the shell number m(t), our above restriction 



(59) 



(60) 



reads 



m(t) <C n x , n 2 . (61) 

We thus may assume that the first two factors in (59) are independent of the last 
and also independent of each other due to eq. (60). On the other hand, the last two 
factors are assumed to be fully correlated, as r ~ r m . With eq. (50) we obtain 



C pini , P2 n 2 ( T ) ~ 



dm 



u ni {t) 


"X 


Un 2 (t) 


U m {t) 




u m (t) 



P2 * 



d/3k f j^- {pi+P2 ^6 [(1 - P(t)) In k m + In • 



(62) 



The shells n\ and m are assumed to be independent. So we only have to plug in the 
definition of the scaling exponents ((p\) } obtaining 



u ni {t) 


Pi \ 


h 


U m {t) 


r 


h 



-C(pi) 



(63) 



Doing the same for the second factor, and performing the integral over m (i.e., 
replacing k m by r -1 ^ 1- ^), we find 

C Pini , P2 n 2 (r) ~ k-^k-^ I d p T -(fW)-( P ^H( Pl )H(P2))/(i-P)_ (g4) 

A saddle point approximation leads to, 

C Pini , P2 n 2 (r) ~ k-^k-^r-^^ (65) 

with 



<t>(pi,P2) = maxp 



f(P)-(pi+P2)P + ((pi) + ((p2) 

1-/3 



(66) 



which is our final result of this section. We think it is worth while to numerically 
check our prediction (66) within the GOY dynamics and possibly try to extend it 
to Navier-Stokes dynamics. 
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4.3 Lyapunov indices and correlation functions 

As one application of the correlation function analysis of the previous subsection, 
we turn to an analysis of the largest Lyapunov index of the GOY model, called Aj,. 
Our work follows closely the approach of A. Crisanti, M.H. Jensen, G. Paladin, and 
A. Vulpiani [49], hereafter called CJPV. 

By defining 8u the infinitesimal difference between two (vector) fields evolving 
under the same equation (i.e. the GOY model in our case), we can define the 
exponential divergence rate 7 after a time delay r of two trajectories close at time 

*' m i 7 IIM* + t)II , fi7 , 

7r (t) = -In . 67 

T \\f>u(t)\\ 

Using (67) we can define Aj, as the limit for infinite r of 7 T (t) and 70 (t) or simply 
7 (t) as the local divergence rate for trajectories. 

CJPV assumed that the Lyapunov exponent should be proportional to the av- 
erage of the largest characteristic rate of change of the system, namely the inverse 
of the eddy turnover time at the dissipation scale [49]. In the GOY model, this 
rate is 

7 = u n k n ~ kn~ a - (68) 

We define n to correspond to a dissipation scale by equating the viscous term D n ~ 
vk 2 n ll n and the nonlinear term C n ~ k n U 2 of the GOY equation, giving 

U n ( t )/k n (t) = v ~ Re' 1 . (69) 

Here we have called v the inverse Reynolds number Re of the GOY dynamical 
equations. From (68) and (69) we obtain 

\ L =< 1 (t) >= J k n 1 - a+na) 6[{l + a)(lnifc„) - \nRe]dadn ~ Re\ (70) 

where 

f(a) + I - a 

a = max„ . I (11 

1 + a 

The numerical value of a estimated by CJPV is 0.46. With f(a) from (56) we 
obtain a = 0.47. K41 predicts a = 0.5. It follows that the exponent a is not 
strongly effected by the intermittency of the system. 

The next step is to calculate the variance in the Lyapunov indices, which CJPV 
estimated by numerically evaluating the integral 

Vl- (ho{t + s)- X L ][ 7o (t)- X L ])ds (72) 
Jo 

which measures the strength of fluctuations in the system. In particular it has been 
shown [36] that ^ = 1 separates weak from strong intermittency. The numerical 
computations [49] show that //£, ~ Re w where w = 0.8. 
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To estimate the integral in equation(72) one notices that the integrand involves 
a correlation of two factors at different times. This is just the kind of integral 

estimated in the previous subsection, except that n is fixed to lie at the dissipation 
scale. Instead of redoing the same calculations it is instructive to use an order of 
magnitude argument to estimate 

When the time separation in equation(72) is of order of the integral scale- value, 
s ~ 1, the fluctuations in *y(t) are also of relative order unity. This result is a 
consequence of the fact that *y(t) is proportional to a velocity, and the velocity is a 
product of roughly independent random variables, one for each shell. Each variable 
has fluctuations of order unity, and the ones for the first shell in *y(t) and *y(t + s) 
are only weakly correlated when s is of order unity. Thus the fluctuations in the 
integrand are of order of the uncorrelated part, so that this part of the integral gives 



where is the dissipative time scale. For s < 1 there is of course no contribution. 
To estimate the remaining part J Td of the integral one goes back to equation(65) and 
imagine integrating a result like this over r. If <f> is greater than one, the integral 
will contribute for small t; otherwise, the main contribution will occur for r of order 
unity. But, using the formula following from She and Leveque, equation (56), one 
can see that <f> is far smaller than one. Hence the main contributions comes at the 
integral scale and we can argue 



For the standard GOY model parameters e = 0.5, A = 2 this gives us an estimate 
fiL ~ Re 2 ' 0A7 = Re ' 94 which has to be compared with the numerical result of C JPV, 
fiL ~ Re ' 8 . Considering the numerical uncertainty and the lack of rigor of our order 
of magnitude argument, the agreement is not too bad. 

We also have to remind that the scaling properties of Xl are linked to small 
fluctuations of the multiplicative process, i.e. in the region where the SL formula 
disagrees with our numerical findings. 

The computations discussed so far make the (strong) assumption that the in- 
stantaneous Lyapunov exponent is controlled by the instantaneous velocity field at 
the dissipation time scale. In fact, between the two quantities there might be a time 
lag tr which characterizes the relaxation time scale for the instantaneous Lyapunov 
exponent to become proportional to the inverse dissipation time scale, tr is not sim- 
ply linked to the scaling properties of the GOY model and its multifractal behavior. 
If tr is small compared to the dissipation time scale, then the analysis previously 
described gives the expected scaling of the intermittent exponent // with respect to 
the Reynolds number. On the other hand, for large enough tr one should find a 
smaller intermittency effect in the system. 

These examples give us some feeling that we might have a crude but useful 
understanding of time dependence in the GOY model. Yet more numerical and 
analytical analyses is definitely necessary. 




(73) 



(74) 
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5 Conclusions 



The dynamical equations of the GOY model permit a solution of a type described 
by BBP in which there are roughly independent fluctuations in all of the shells 
of the inertial range. The long-ranged structure of these fluctuations give all the 
multifractal scaling properties. The conservation laws for energy and helicity play 
an essential role in the structure of the solution, but such 'details' as the momen- 
tum scale A, the Reynolds number, and the form of the viscous cutoff might also 
determine the exact scaling exponents. 

All scaling that we have examined in the model seems to be determined by one 
set of multifractal exponents, e.g. those of the velocity. In particular, many time- 
dependent correlations may be estimated from these exponents. However, at the 
moment there exists no accurate way of calculating these scaling exponents except 
by direct numerical simulation. However, see [41] and [50] for an interesting set of 
insights into the possible structure of the multifractal behavior. 

This paper devotes some attention to the scaling of time correlations. The cor- 
responding frequency spectra have recently been studied in [51]. Notice that these 
are the quantities which are measured in experimental studies of turbulence. Turbu- 
lence theory up to now has mainly focused on spatial structure functions and wave 
vector spectra and has connected them with time structure functions and frequency 
spectra only via Taylor's hypothesis. Clearly time-dependence is worth studying in 
its own right. 

We finally stress, that the exact way multifractality works itself out in the GOY 
dynamics is slightly different from what happens in real turbulence because there 
is no analog of the sweeping which plays such an essential role in Navier- Stokes 
dynamics, and, of course, because the GOY equations are only a model which might 
or might not catch essential features of the Navier-Stokes equations. 
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Tables 



n 


numerics 


eq (25) 


numerics 


eq (25) 




e = 0.3 


e = 0.3 


e = 0.5 


e = 0.5 


4 


15000 


15000 


0.14 


0.14 


5 


1.31664 


1.31667 


1.35 


1.10 


6 


0.28291 


0.28291 


0.27 


0.36 


7 


0.76848 


0.76857 


0.78 


0.59 


8 


0.37766 


0.37770 


0.37 


0.46 


9 


0.61321 


0.61333 


0.62 


0.52 


10 


0.43483 


0.43532 


0.43 


0.48 


11 


0.55055 


0.55110 


0.55 


0.50 



Table 1: The r n from our numerical calculation, compared with the expression (25) 
when e = 0.3 and e = 0.5. In the latter case, we show the time averaged veloci- 
ties. Note that the same alternation between high and low values occurs in all four 
columns. 



q 


derived from ( q via (55) 


direct numerical result 


derived from (32) 


1.0 


0.00 


0.0005 ± 0.0002 


0.00 


1.5 


0.41 


0.42 ±0.03 


0.40 


2.0 


1.00 


1.00 ±0.03 


1.00 


2.5 


1.72 


1.71 ±0.04 


1.77 


3.0 


2.54 


2.52 ±0.1 


2.67 


3.5 


3.44 


3.41 ±0.2 


3.68 


4.0 


4.42 


4.25 ±0.2 


4.78 


4.5 


5.44 


5.23 ±0.4 


5.95 



Table 2: The ratio fi q / fi 2 from our data for velocity scaling compared with the ratio 
deduced from the the scaling analysis of dissipation. The numerical errors are esti- 
mated from the differences between the fits in slightly different fitting regions. The 
error in the first column due to the numerical error in the ( q is much less significant 
compared to the error in the second column. In the last column we give fi q / fi 2 as it 
follows from eqs. (32) and (45). These values are systematically slightly larger than 
the first two columns, yet the agreement is satisfactory. 
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Figure captions 



Figure 1: Elimination of period three oscillations. We show three curves plotted 
against shell number n. Two are for s nj i = (\U n \): A plot for the normal viscosity 
case (D n = vkfjJn with cf> = 2) and another for the hyperviscosity case ((f) = 4). 
Both of these show period-three oscillations. The remaining curve plots S n> i for 
the hyperviscosity case and shows no period-three oscillations. In each case time 
averages are taken over 15 000 large eddy turnover times, the number of shells is 
22. One curve is for our standard model parameters, the hyperviscosity cases have 
v = 10" 15 . 



Figure 2: Different kinds of scaling analyses with increasing accuracy: s, which is 
the magnitude of the velocity, S which is a cube root of the imaginary part of a 
product of three velocities, and E, which is the cube root of the energy flux. Part a 
shows curves drawn for q = 1 ; part b for q = 6. In both cases, E gives the longest 
scaling range, and hence probably the best estimates for scaling exponents. All 
curves are drawn for the standard parameter values. For comparison, K41 scaling 
is also shown. 



Figure 3: The deviations 6( q from Kolmogorov scaling are compared with one an- 
other for different parameter sets. We found it most instructive to plot S( q /(q(q — 3)) 
against q, to eliminate the trivial agreement at q = and q = 3. The dashed line 
is our numerical calculation, with the standard parameter values. The solid line, 
which is impressively close to our numerical results for q > 3, is the result (32) of 
She and Leveque [41]. The dotted line reports the shell results for another set of 
parameters, namely e = 0.7 and A = 2. Note the pronounced disagreement with the 
SL result in this case. 
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Figure 4: 6( q versus q for four sets of parameter values. Three parameter pairs 
(e, A) lie on the curve (33) which defines the right value of the helicity. These have 
A = 10/3 and A = 10/7 paired with their corresponding e's. The last value lies off 
the curve and has e = 0.7 with A = 2. Note how the values on the curve stand 
grouped together in comparison with the other one. 



Figure 5: Scaling corrections 6( q for three different forms of viscous damping. In 
each case, the damping term is D n = vk n ^U n . The solid line shows the result for 
<f) = 2 and v = 10~ 7 , the dashed line for <f> = 4 with v = 10~ 15 , and the dotted line 
for <f> = 6 with v = 5 • 10~ 25 . 



Figure 6: The top picture shows the scaling behavior of e Tj2 in r, while the bottom 
employs the idea of extended-self similarity to plot one moment against another, i.e. 
e Tj4 against e Tj2 . The slope of the best fit is also shown (dashed line). 
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